package tw;

abstract class rk4 {

	public double dydx[];

	public double y[];

	public double yout[];

	int n;

	rk4() {
	}

	public abstract void derivs(double d, double ad[], double ad1[]);

	public void init(int dim, boolean self) {
		n = dim;
		dydx = new double[n];
		yout = new double[n];
		if (self)
			y = yout;
		else
			y = new double[n];
	}

	public void moveNext(double h) {
		core(0.0D, h);
	}

	public void moveNext(double h, double yin[]) {
		y = yin;
		core(0.0D, h);
	}

	private void core(double x, double h) {
		double dym[] = new double[n];
		double dyt[] = new double[n];
		double yt[] = new double[n];
		double hh = h * 0.5D;
		double h6 = h / 6D;
		double xh = x + hh;
		derivs(xh, y, dydx);
		for (int i = 0; i < n; i++)
			yt[i] = y[i] + hh * dydx[i];

		derivs(xh, yt, dyt);
		for (int i = 0; i < n; i++)
			yt[i] = y[i] + hh * dyt[i];

		derivs(xh, yt, dym);
		for (int i = 0; i < n; i++) {
			yt[i] = y[i] + h * dym[i];
			dym[i] += dyt[i];
		}

		derivs(x + h, yt, dyt);
		for (int i = 0; i < n; i++)
			yout[i] = y[i] + h6 * (dydx[i] + dyt[i] + 2D * dym[i]);

	}
}
